arXiv: 1503.01109v2 [astro-ph.CO] 9Jul2015 


Mon. Not. R. Astron. Soc. 000, 1-14 (2015) Printed 10 July 2015 (MN BTgX style file v2.2) 


Exploring the liminality: properties of haloes and subhaloes in 
borderline f(R) gravity 


Difu Shi 1 *, Baojiu Li 1 , Jiaxin Han 1 , Liang Gao 2,1 , Wojciech A. Hellwing 1 

1 Institute for Computational Cosmology, Department of Physics, Durham University, South Road, Durham DHI 3LE, UK 

2 National Astronomical Observatories, Chinese Academy of Sciences, 20A Datun Road, Chaoyang District, Beijing 100012, China 


10 July 2015 


ABSTRACT 

We investigate the properties of dark matter haloes and subhaloes in an f(Jl) gravity model 
with l/ijol = 10 using a very high-resolution N-body simulation. The model is a border¬ 
line between being cosmologically interesting and yet still consistent with current data. We 
find that the halo mass function in this model has a maximum 20% enhancement compared 
with the ACDM predictions between z = 1 and z = 0. Because of the chameleon mechanism 
which screens the deviation from standard gravity in dense environments, haloes more mas¬ 
sive than 10 13 /i _ 1 Mq in this f(R) model have very similar properties to haloes of similar 
mass in ACDM, while less massive haloes, such as that of the Milky Way, can have steeper 
inner density profiles and higher velocity dispersions due to their weaker screening. The halo 
concentration is remarkably enhanced for low-mass haloes in this model due to a deepening 
of the total gravitational potential. Contrary to the naive expectation, the halo formation time 
Zf is later for low-mass haloes in this model, a consequence of these haloes growing faster 
than their counterparts in ACDM at late times and the definition of Zf. Subhaloes, especially 
those less massive than 10 11 h~ 1 Mo, are substantially more abundant in this f(R) model for 
host haloes less massive than l0 13 h~ 1 M@. We discuss the implications of these results for 
the Milky Way satellite abundance problem. Although the overall halo and subhalo properties 
in this borderline f(R) model are close to their ACDM predictions, our results suggest that 
studies of the Local Group and astrophysical systems, aided by high-resolution simulations, 
can be valuable for further tests of it. 

Key words: gravitation - methods: numerical - galaxies: halos - cosmology: theory - dark 
matter - large-scale structure of Universe 


1 INTRODUCTION 

The observed accelerated cosmic expansion is one of the most puz¬ 
zling problems in modern physics (e.g., Weinberg et al. 2013). In 
less than twenty years, it has motivated the proposal of a huge num¬ 
ber of models. Apart from the current standard ACDM paradigm, in 
which the acceleration is driven by a cosmological constant A, such 
models are divided roughly in two classes. The first class introduces 
new physics in the particle sector and suggests that the acceleration 
is due to some new matter species, often known as dark energy. The 
second class proposes new physics in the gravity sector, so that the 
standard theory of gravity, Einstein’s General Relativity (GR), is 
modified on cosmological scales to accommodate the accelerated 
expansion. This latter class of theories are commonly referred to 
as modified gravity (Clifton et al. 2013; Joyce et al. 2014) which is 
increasingly becoming an active research area. 

For over a decade, the well-known f(R) gravity (Carroll et al. 
2004, 2005) model has been a leading modified gravity candidate 
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to explain the cosmic acceleration, although it actually has a much 
longer history in other contexts. It is a subclass of the more general 
theory called the chameleon theory (Khoury & Weltman 2004), in 
which an extra scalar degree of freedom is invoked, that can medi¬ 
ate a modification to the standard gravitational force of GR (known 
as the fifth force). This deviation from GR is not necessarily ruled 
out by local experiments, as the theory can employ the chameleon 
mechanism (Khoury & Weltman 2004) to suppress the fifth force in 
dense environments such as the Solar System (see below for more 
details about this so-called chameleon screening). This means that 
the theory could pass local gravity tests. However, in less dense en¬ 
vironments, such as those encountered on cosmological scales, the 
deviation from GR becomes sizeable, which means that cosmology 
can provide a unique means to probe new physics of this kind. 

There are several important features of f(R) gravity, some of 
which seem to have not been emphasised enough. Firstly, it is a well 
accepted perception that f(R) gravity is flexible and, thanks to its 
4th-order field equations, can, in principle, accommodate arbitrary 
background cosmologies (see, for instance, He & Wang 2013, for 
a concrete example of ACDM background cosmology). In spite of 
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the general impression that ‘ f(R ) gravity can accelerate the cos¬ 
mic expansion’, it should be noticed that there is no necessary con¬ 
nection between the ‘acceleration’ and ‘modified gravity’ parts of 
f(R) gravity: the well-studied model of Hu & Sawicki (2007), as 
an example, can essentially be written as a cosmological constant 
plus a modification to the GR gravitational law. In this sense, f(R) 
gravity is not a ‘better’ model than ACDM, but rather one with a 
different nature of gravity, the study of which can shed light on the 
question why GR is successful and whether cosmological data can 
disapprove it. 

Secondly, the chameleon screening in f(R) gravity is a mech¬ 
anism that ‘could work’, but not necessarily ‘will work’ - whether 
it works depends on the system under consideration and the func¬ 
tional form of f(R ). Again, taking the Hu & Sawicki (2007) model 
for example: how efficiently the screening works is determined by 
a model parameter | fno | (and another parameter n which is often 
fixed; see below). Increasing this parameter makes it less likely for 
the screening to be effective. Values of |/ho| < 1CD 6 are more dif¬ 
ficult to distinguish from GR using cosmological observations, im¬ 
plying a limit of cosmological constraints. On the other hand, there 
are recent claims that |/ho| = 10~ 6 could be in tension with as- 
trophysical observations (see, e.g., Jain et al. 2013). Given that the 
chameleon mechanism works with different efficiency in different 
environments, it is critical to examine whether these stringent con¬ 
straints become weaker when the environments of the astrophysi- 
cal systems are more accurately modelled. The same could be said 
about terrestrial tests of the chameleon theory (see, e.g., Brax et al. 
2007a,b, for references to some pioneering works in this direction). 

For this reason, the Hu & Sawicki (2007) f(R) gravity model 
with |/flo| = 10 -6 could be considered as being borderline be¬ 
tween cosmological and astrophysical constraints: for higher values 
the model will probably have trouble with local and astrophysical 
tests, and for lower values the model is likely to be no longer inter¬ 
esting cosmologically. Here, we suggest a ‘bisection’ approach to 
the study of f(R) gravity: we first conduct a detailed investigation 
of the cosmological and astrophysical implications of the model 
with |/ro| = 1CP 6 , and then push the study and the resulting con¬ 
straints to larger or smaller values based on the outcome. We hope 
to use this paper, in which we will concentrate on the cosmological 
aspects, as an initial step in this direction, to motivate further, more 
in-depth, studies. 

In this work, we employ one of the highest-resolution N-body 
simulations of f(R) gravity currently available to study its effects 
on the properties of dark matter haloes and their subhaloes. These 
are the fundamental building blocks of the large-scale structure of 
our Universe and are closely connected with cosmological observa¬ 
tions such as galaxy surveys. Previous studies have shown that the 
model considered here makes rather similar predictions to GR for 
many other cosmological observables, such as the matter and veloc¬ 
ity power spectra (Li et al. 2013; Hellwing et al. 2013; Zhao 2014; 
Taruyaetal. 2014), void properties (Cai et al. 2014; Zivicketal. 
2014), redshift space distortions (Jennings et al. 2012), the inte¬ 
grated Sachs Wolfe effects (Cai et al. 2014) and X-ray scaling rela¬ 
tions of clusters (Arnold et al. 2014a), but the simulation resolution 
used have not been high enough to study haloes and subhaloes in 
great detail (see, e.g., Corbett Moran et al. 2015, for a recent high- 
resolution zoom-in simulation which has a different focus from that 
of this paper). 

This paper is structured as following: In §2 we very briefly 
describe the f(R) model studied here and summarise the technical 
specifications of our simulations. §3 and §4 present our detailed 
analyses of halo and subhalo properties respectively, and compar¬ 


isons with the ACDM model. Finally, we summarise and conclude 
in §5. 

Throughout this paper, we use the unit c = 1, where c is the 
speed of light. 


2 f(R) GRAVITY AND SIMULATIONS 

In this section we briefly review the general theory of f(R) gravity 
(§2.1), motivate the model which we focus on (§2.2) and describe 
the algorithm and technical specifications of our cosmological sim¬ 
ulations (§2.3). 


2.1 f(R) gravity and chameleon screening 


The f(R) gravity model is designed as an alternative to dark energy 
to explain the accelerated expansion of the Universe. It generalises 
the Ricci scalar R to a function of R in the Einstein-Hilbert action, 

s = <» 

where G is Newton’s constant and g is the determinant of the metric 

5 m v 

Minimising the action Eq. (1) with respect to the metric tensor 
leads to the modified Einstein equation 


G ^iir + f rR/xii g fit/ 


2 f ~ nf * 


V^„f(R) = 8nGT™, ( 2 ) 


where G /JU is the Einstein tensor, fn = df /dR, V M is the covari¬ 
ant derivative, □ = V“V Q and T^l the energy momentum tensor 
for matter fields. As R contains second-order derivatives of g^, 
Eq. (2) has up to fourth-order derivatives. It is helpful to consider it 
as the standard Einstein equation for general relativity with an ad¬ 
ditional scalar field fn. By taking the trace of Eq. (2), the equation 
of motion for frt can be obtained as 

n/i? = — (R — /rR + 2/ + 8nGpm ), (3) 

where p m is matter density. 

We consider a flat universe and focus on scales well below the 
horizon. On these scales, we can apply the quasi-static approxima¬ 
tion by neglecting the time derivatives of /r? in all field equations 
(see, e.g., Bose et al. 2015, for tests which show that this approxi¬ 
mation works well for the model studied here). Then Eq. (3) sim¬ 
plifies to 

V 2 /fl = [-R(/h) - R + 8nG(pm - pm)] , (4) 

in which V is the three-dimensional gradient operator and an over¬ 
bar means we take the cosmological background value of a quan¬ 
tity. a is the cosmic scale factor, normalised to a = 1 today. Simi¬ 
larly, the modified Poisson equation, which governs the Newtonian 
potential 4? in f(R) gravity, can be simplified to 

V 2 $ = ±^a 2 ( Pm - p m ) + i [R{f R ) - A] • (5) 

There are two distinct regimes of solutions to the above equa¬ 
tions: 

• when \fn\ <g |4>|. the GR solution R = —8nGp m holds 
to a good approximation and one has V 2< f> « AivGSpm where we 
have defined Sp m = p m — p m , as the matter density perturbation. 
The effect of modified gravity is suppressed in this regime, which is 
a consequence of the scalar field being screened by the chameleon 
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mechanism (Khoury & Weltman 2004). 

• when \fu\ ^ |$|, one has |5i?| <C <5p m where SR = R — R, 
and so V 2< f> ss 16nG8p m /S. Compared with the standard Poisson 
equation in GR, we see a 1/3 enhancement in the strength of grav¬ 
ity regardless of the functional form of f(R). This is known as the 
unscreened regime, in which the chameleon mechanism does not 
work efficiently. 

The chameleon mechanism is so named because it is most ef¬ 
ficient in dense environments (or, more precisely speaking, regions 
of deep gravitational potential), where the scalar field fn acquires 
a heavy mass and the (Yukawa-type) modified gravitational force it 
mediates decays exponentially with distance so that it cannot be de¬ 
tected experimentally. The Solar System is one example of such an 
environment where f(R) gravity might be in the screened regime 
and thus viable (i.e., not yet ruled out by local gravity experiment). 
However, to determine whether a specific f(R) model is indeed vi¬ 
able is much more difficult, because this depends on the large-scale 
environments of the Solar System, such as the Milky Way Galaxy 
and its host dark matter halo. To assess this therefore requires high- 
resolution numerical simulations that can accurately describe these 
environments, and this is one goal of our paper. On the other hand, 
even if an f(R) model passes local tests, there is still a possibility 
that it deviates significantly front GR on cosmic scales, where the 
chameleon mechanism is not as efficient. To study such deviations 
also requires accurate numerical simulations. 


2.2 The f{R) model of this work 


In this work we study the model proposed by Hu & Sawicki (2007, 
hereafter HS), which is specified by the following functional form 

of f(R): 


f(R) = ~M 2 


ci ( -R/M 2 ) n 
c 2 {-R/M 2 ) n + 1 ’ 


(6) 


in which ci, C 2 are dimensionless model parameters, and M 2 = 
8nGp m o/3 = HoQ, m is another parameter of mass dimension 2; 
here H is the Hubble rate and H m is the present-day matter en¬ 
ergy density in units of the critical density (p c = 3H 2 /8nG). We 
always use a subscript o to denote the current value of a quantity 
unless otherwise stated. 

When the background value of the Ricci scalar satisfies |ii| 

M 2 , we can simplify the trace of the modified Einstein equation of 
this model as 


— R 


8tt Gpm - 2/ 



(V) 


This is approximately what we have for the background cosmology 
in the standard ACDM model, with the following mapping 


ci Ha 

— = 6 7S — > 

C2 * lm 


( 8 ) 


where Ha = 1 — fi m . 

By taking Ha « 0.7 and fi m ss 0.3, we have \R\ ss 40 M 2 
M 2 today (remember that |.R| is even larger at earlier times), and 
so the above approximation works well. Moreover, this can be used 
to further simplify the expression for fn\ 


fR 




< 0. 


(9) 


This can be easily inverted to obtain R(fn), which appears in the 
scalar field and modified Poisson equations as shown above. As a 
result, two combinations of the three HS model parameters, namely 


Table 1. The parameters and technical specifications of the N-body simula¬ 
tions of this work. e s is the threshold value of the residual (see, e.g., Li et al. 
2013, for a more detailed discussion) for the convergence of the scalar field 
solver. Note that N Te f is an array because we take different values at differ¬ 
ent refinement levels, and that erg is for the ACDM model and only used to 
generate the initial conditions - its value for f(R) gravity is different but is 
irrelevant here. 


parameter 

physical meaning 

value 


present fractional matter density 

0.281 

Ha 

1 

0.719 

Hf, 

present fractional baryon density 

0.046 

h 

#o/(100 km s -1 Mpc -1 ) 

0.697 

n s 

primordial power spectral index 

0.971 


r.m.s. linear density fluctuation 

0.820 

n 

HS f(R) parameter 

1.0 

Sro 

HS f(R) parameter 

-1.0 x 10“ 6 

-^box 

simulation box size 

64 h _1 Mpc 

N p 

simulation particle number 

512 3 

m p 

simulation particle mass 

1.52 x 10 8 h~ 1 M Q 

N dc 

domain grid cell number 

512 3 

Aref 

refinement criterion 

3, 3,3, 3,4, 4, 4,4... 

e s 

scalar solver convergence criterion 

10" 8 

ef 

force resolution 

1.95 L -1 kpc 

-^snap 

number of output snapshots 

122 

+ni 

redshift when simulation starts 

49.0 

■^final 

redshift when simulation finishes 

0.0 


n and ci/c§, completely specify the model. In the literature, how¬ 
ever, this model is often specified by fuQ instead of ci/cl, because 
fno has a more physical meaning (the value of the scalar field to¬ 
day), and the two are related by 


Cl 


1 

n 


3 1 + 4 


Ha 

Hm 


n+1 

fHo¬ 


rn 


We will focus on a particular HS f(R) model with n = 1 and 
I/ho | = 1(T 6 , which is sometimes also referred to in the literature 
as F6. Such a choice of /no is made deliberately as a borderline : 
models with |/no| + ICR 5 are likely to already be in tension with 
cosmological observations (see, e.g., Lombriser 2014, for a review 
of current constraints on f(R) gravity), while those with |/no| < 
10 -6 are generally hard to distinguish from ACDM. 


2.3 Cosmological simulations of f(R) gravity 

The simulation used in this work was executed using the ECOSMOG 
code (Li et al. 2012). ECOSMOG is a modification to the publicly 
available N-body and hydro code RAMSES (Teyssier 2002). New 
routines were added to solve the scalar field and modified Einstein 
equations in f(R) gravity. This is a massively parallelised adaptive 
mesh refinement (AMR) code, which starts off from a uniform grid 
(the so-called domain grid) covering the cubic simulation box with 
N d ' c cells on each side. When the effective particle number in a 
grid cell exceeds a pre-defined criterion (iV re f), the cell is split into 
8 daughter cells so that the code hierarchically achieves higher res¬ 
olutions in dense environments. Such high resolutions are needed 
both to accurately trace the motion of particles and to ensure the ac¬ 
curacy of the fifth force solutions. The force resolution, et, is twice 
the size of the cell which a particle is in, and we only quote the 
force resolution on the highest refinement level. 
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Figure 1. Comparison of differential halo mass functions in GR (red circles) and F6 (blue triangles) with the Sheth & Tormen (1999) prediction for GR, at 
three redshifts — z = 0.0 (left panel), 0.5 (middle panel) and 1.0 (right panel). The relative difference between the two models is plotted in the bottom panels. 
Haloes are identified using a FoF algorithm with linking length 0.2. As we only have one realisation, the error bars are estimated from subsampling by dividing 
the simulation box into eight subboxes of equal size; the difference between F6 and GR mass functions, Ani, was computed for each subbox i, and its mean 
value ((An)) and standard deviations (cresn) were obtained using the values from the 8 subboxes; the relative difference was then calculated as (An)/(n^i; ), 
with the error bars obtained in the standard way of error propagation. The vertical dashed line indicates a cut of our FoF halo catalogue at ~ 700 particles, or 
MpoF ~ 1 0 11 h ~ 1 Mq , for illustrative purpose, above which there is good agreement between the GR mass functions and the Sheth-Tormen fitting formulae 
(better than 10%). 


The parameters and technical specifications of our simulations 
are listed in Table 1. The cosmological parameters are adopted from 
the best-fit ACDM cosmology of WMAP9 (Hinshaw et al. 2013). 
The simulation was evolved from an initial redshift Zi n i = 49 to to¬ 
day, and the initial conditions were generated using the MPGRAFIC 
package (Prunet et al. 2008). For comparison, we ran a f(R) and a 
ACDM simulation using exactly the same initial conditions and the 
same technical specifications (we have used the ACDM initial con¬ 
ditions for the f(R) gravity simulation because these two models 
are practically indistinguishable at epochs as early as 2 i n i = 49). 
The small size of our simulation box implies that the properties of 
high-mass objects, such as their number densities, could be subject 
to run-by-run variations. However, the fact that our f(R) and GR 
simulations start from the same initial conditions helps to suppress 
the run-by-run variation when we look at the relative difference be¬ 
tween the predictions of the two models. 


With 512 3 particles in a box of size Lbox = 64Ii _1 Mpc, this 
is currently the highest resolution cosmological simulation of f(R) 
gravity which runs to z = 0. Another high-resolution f(R) simula¬ 
tion has been conducted (Corbett Moran et al. 2015), in which the 
zoom-in technique was used to study the effects of f(R) gravity on 
a Virgo-cluster-scale dark matter halo. Both simulations are purely 
dark matter. Recently, a hydrodynamical simulation was carried out 
by Arnold et al. (2014b), which had a higher particle resolution and 
focused mainly on a different model parameter and early times, at 
which the model studied here is almost indistinguishable from GR. 


3 PROPERTIES OF DARK MATTER HALOS 

In this section we will concentrate on the properties of dark matter 
haloes measured from our simulations. Dark mater haloes are the 
most basic blocks of the large-scale structure and host the formation 
and evolution of galaxies. Therefore, the study of their properties is 
of great importance to the understanding of the fundamental nature 
of gravity. A number of halo properties have been studied in de¬ 
tail in the context of f(R) gravity, such as the angular momentum, 
spin, velocity dispersion (Lee et al. 2013; He et al. 2015), veloc¬ 
ity profile (Gronke et al. 2014) and screening (Zhao et al. 2011a; 
Li et al. 2012; He et al. 2014). The improved resolution of our sim¬ 
ulations enables us to study a wider range of the physical properties 
of haloes. 

3.1 Halo mass functions 

The differential mass function, dn(M, z)/d log M, defined as the 
number of dark matter haloes per unit logarithmic mass found per 
unit volume, is an important theoretical and observational statistic 
of the dark matter density field. Indeed, the abundance of dark mat¬ 
ter haloes is sensitive to the underlying cosmological model. Both 
N-body simulations and (semi-)analytical formulae have been used 
to predict the halo mass function (see, e.g., Sheth & Tormen 1999; 
Jenkins et al. 2001; Reed et al. 2007, for some examples of analyt¬ 
ical mass function fitting formulae). 

In order to compare with the above-mentioned fitting formu¬ 
lae, we use the friends-of-friends (FoF) group-finding algorithm to 
identify dark matter haloes, using a linking length of 0.2 times the 
mean inter-particle separation (Davis et al. 1985). 

In Fig. 1, we plot the differential halo mass function measured 
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from our simulations, along with the theoretical prediction for GR 
from Sheth & Tormen (1999) (upper panels), and the relative differ¬ 
ence between f(R) gravity and GR (lower panels) at z = 0 (left), 
0.5 (middle) and 1.0 (right). For the mass range we consider, the 
Sheth & Tormen (1999); Jenkins et al. (2001) & Reed et al. (2007) 
fitting formulae all agree reasonably well and so we only plot one 
of them. We can see from the upper panels that the fitting formula 
describes very well the FOF halo mass function for GR at the red- 
shifts studied, down to a halo mass of about 2 — 3 x 10 10 Ii~ 1 Mq 
(which corresponds to ~200 simulation particles). The mismatch 
at masses above ~ 10 14 /i -1 Mq is due to the lack of volume for 
our small simulation box. 

Fig. 1 (lower panels) indicates that the differential halo mass 
function for F6 model studied here is up to ~20% larger than the 
result for a ACDM model with the same cosmological parameters. 
The difference is purely a result of the modified gravitational force 
in the F6 model. However, due to the strong chameleon screening in 
this model, the enhancement is very mild and hard to detect obser- 
vationally. This is why we call F6 a borderline model - it probably 
represents the limit achievable by many cosmological observations 
for the near future, even though it might still potentially be ruled out 
by employing certain observables (e.g., Schmidt 2010; Zhao et al. 
201 lb; Bel et al. 2014; Lombriser et al. 2015), or using astrophys- 
ical observations (e.g. Jain et al. 2013). 

Inspecting the lower panels of Fig. 1 more closely, we observe 
the trend that for very massive haloes, the mass functions for f(R) 
gravity and GR agree, which is because the chameleon mechanism 
works efficiently for such haloes to suppress the effects of modi¬ 
fied gravity. Disagreement between the two models appears below 
some critical mass, which increases with time, because at late times 
the chameleon mechanism is less efficient at suppressing modified 
gravity. Finally, at very low halo masses, we see the trend that GR 
starts to produce more haloes than /(/?.) gravity, which is a result 
of a larger fraction of small haloes having been absorbed into big 
haloes in f(R) gravity (Li & Efstathiou 2012). 

It is well known that certain properties of dark matter haloes, 
such as the mass function, depend on the halo definition used (e.g.. 
White 2001; Sawala et al. 2013). In the above, to make comparison 
with the Sheth-Tormen formulae, we have used FOF haloes. When 
studying halo properties, what is more often used in the literature 
is M200. the mass inside the radius r2oo within which the average 
density is 200 times the critical density, p c . To check whether the 
choice of the halo definition affects our result, we plot in Fig. 2 the 
difference between the f(R) and GR mass functions when using 
A/ 200 , again at z = 0.0 (upper panel), 0.5 (middle) and 1.0 (lower 
panel). We find the same qualitative features as in the lower pan¬ 
els of Fig. 1, but also some quantitative differences in the curves. 
In particular, the curves are smoother and better-behaved when us¬ 
ing M 200 , which may be because the FOF haloes are too irregular 
in their shapes and gravity is enhanced with different efficiency in 
different parts of the haloes, which can contaminate the screening 
effect expected for ideal spherical haloes (see, e.g., Li & Efstathiou 
2012; Li & Lam 2012; Lombriser et al. 2013, 2014, for more dis¬ 
cussion about the expected behaviour of the f(R) halo mass func¬ 
tion). 

We will use M 200 in the rest of this paper, because of its wide 
use in the literature. Furthermore, to ensure good resolution of halo 
structure, we will conservatively restrict our analysis to haloes with 
more than 700 particles (M200 > 10 11 h^ 1 M@). Even cut at ~400 
particles, we have found that the FoF mass functions at z = 0, 0.5 
and 1 show agreements with fitting formulae better than 10%. 
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Figure 2. The ratios of the differential halo mass functions between /(fij 
gravity and GR, for the same three redshifts as in Fig. 1. Here the halo mass 
is M 200 , defined as the mass within the radius at which the average density 
200 times the critical density. The error bars are calculated in the same 
way as in Fig. 1. The vertical dashed line indicates roughly the smallest 
halo mass (700 particles, or M 200 ~ 10 11 h ~ 1 Mq) we have used in the 
analyses of this paper. 


3.2 Mass distribution inside haloes 

The inner structure of dark matter haloes provides invaluable infor¬ 
mation about their formation history, which can also be affected by 
the nature of gravity. In this subsection, we look at the dark matter 
density profiles and concentration-mass relations for haloes in the 
two models. 

In Fig. 3, we show the stacked halo dark matter density pro¬ 
files at three redshifts z = 0.0 (left), 0.5 (middle) and 1.0 (right). 
The distances from halo centres, as plotted on the horizontal axis, 
are rescaled by r2oo, and all haloes with M ^ 10 11 hr 1 Mq in our 
simulations are divided into 5 mass ranges as indicated in the leg¬ 
end (the highest mass bin does not show up in the z = 1.0 panel, 
since at that time the very massive haloes have not formed in great 
numbers). Note that the widths of mass bins are different, and we 
do not make finer subdivisions of the three most massive bins since 
the model differences are small there, an observation we discuss 
now. 

From Fig. 3, we see that the density profiles of haloes more 
massive than 10 13 Ii~ 1 Mq show almost no difference between the 
two models at all three redshifts, because these haloes are very 
efficiently screened by the chameleon mechanism. The haloes in 
the mass bin 10 11 ~ 10 12 Ii~ 1 Mq have up to 60% higher den¬ 
sity towards their centres in f(R) gravity than in GR, because 
the screening efficiency is weaker. Thus, Milky-Way-sized haloes 
have steeper inner profiles in f(R) gravity. Note, however, that 
as the force resolution of our simulations is ~ 2 /i _1 kpc, we 
show the results within 5 times of it, i.e., ~ 10 hr 1 kpc, us¬ 
ing open rather than filled symbols. We have explicitly checked 
that 10fi -1 kpc is roughly equal to the Power convergence radius 
(Power et al. 2003; Schaller et al. 2014) in our smallest halo mass 
bins (10 11 ~ 10 12 /i _1 Mq), and is larger than the convergence ra¬ 
dius for other halo mass bins shown in Fig. 3. Though the Power ra¬ 
dius is found by testing convergence on simulations with tree (and 
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Figure 3. Stacked dark matter halo density profiles for five mass bins (indicated by the legend) at three epochs: 2 : = 0.0 (left panel), 0.5 (middle panel) and 
1.0 (right panel). Upper panels : density profiles - results from the GR simulation are shown as circles while the f(R) results are shown as triangles; the solid 
curves connect the symbols and the halo mass increases from the bottom curves to the top curves. Lower panels : the corresponding ratios between f(R) gravity 
and GR which deviate from unity more for smaller mass bins. Open symbols are used when the distance from halo centre is smaller than the force resolution 
(which happens only in the lowest mass bin). The error bars are 1-cr standard deviations for all haloes in each radius bin. To assist visualisation, in the top 
panels, we have rescaled the stacked density profiles of haloes within the mass bins 10 14 ~ 10 15 /i — 1 M@, 10 12 ~ 10 13 h —1 Mq, 3 x 10 11 10 12 /i _1 Mq 

and 10 11 ~ 3 x 10 11 h~^M q by factors of 10, 0.1, 10 — 2 and 10“ 3 respectively. The numbers of haloes in each bin, starting from the most massive one, 
are respectively (numbers for the F6 simulation are in parentheses): 7 (7), 72 (78), 200 (232), 509 (586), 1558 (1714), 3758 (3936) at z = 0, 3 (3), 62 (62), 
194 (215), 539 (624), 1665 (1968), 4224 (4554) at z = 0.5, and 0 (0), 48 (49), 152 (156), 508 (564), 1730 (2047), 4404 (5177) at z = 1.0. 
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Figure 4. Dark matter halo concentrations C200 as a functions of M200 at three redshifts, z = 0.0 (left panel), 0.5 (middle panel) and 1.0 (right panel). 
Results for GR (F6) are shown using red circles (blue triangles), and the curves are power-law fits of the C-M200 relations: in the case of GR (red solid line), 
the relation can be fitted using a single power law for the whole mass range, while for F6 this is no longer true and so we do not show any fitting. The error 
bars are obtained as the 1-cr standard deviation of all haloes in each given mass bin. 
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not AMR) codes, the physics of collisional relaxation used in its 
derivation is the same in our simulations, and so we use it as a ref¬ 
erence. We conclude therefore that the region within 0.1 x r 2 oo of 
Milky-Way-sized (and smaller) haloes in f{R) gravity has a signif¬ 
icantly steeper density profile than in GR. and will further confirm 
this by studying the halo concentration-mass relations below. 

By comparing the three panels in Fig. 3, it is also evident that 
the differences in the inner density profiles of the two models grow 
in time. This is as expected because the effect of modified gravity is 
cumulative, and also because at late times the chameleon screening 
is weaker in general, which leads to stronger modifications to GR. 
In particular, haloes with masses below ~ 3 x 10 11 /i _1 Mq already 
show signifiant discrepancy between F6 and GR at 2 = 1, and for 
haloes front the mass bin 3 x 10 11 ~ 10 12 /i _ 1 M© the discrepancy 
starts at later times because of more efficient chameleon screening, 
although by 2 = 0 the model differences have become roughly the 
same for these two bins. 

Next, we fit the dark matter density profiles in the two models 
using the Navarro-Frenk-White (Navarro et al. 1996, 1997. NFW) 
formula, which is given by 


p{r) 


_Ps_ 

(r/r s )(l + r/r s ) 2 ’ 


( 11 ) 


in which p s and r B are the scale density and scale radius of the halo. 
The p B and r s parameters are connected to the halo mass, AT 200 (or 
equivalently, the virial radius r 200 ), and concentration, c (note that 
we have neglected the subscript in C 200 for brevity), through 


200p c __ 

3 [ln(l + c) — c/(l + c)] ’ 


( 12 ) 


c = r 2 oo/r s . (13) 

In practice, we obtain the M 200 and r 200 of each halo accord¬ 
ing to the spherical over-density definition, and estimate c using 
Eq. (13) from the best-fitting r s . Lombriser et al. (2013) found 
that haloes in f(R ) gravity can be well described by the NFW 
formula Eq. (11). In this work, we have further confirmed this by 
explicitly checking the \ 2 goodness-of-fit, in which we found that 
Eq. (11) works almost equally well in GR and f(R) gravity (with 
marginally smaller \ 2 for haloes between ~ 10 12 — 10 13 h~ 1 Mq 
in f(R) gravity), though the concentration parameters can be dif¬ 
ferent, as we shall show below. 

Fig. 4 shows the halo concentration-mass relation, c(M 2 oo), 
also at three redshifts 2 = 0.0 (left), 0.5 (middle) and 1.0 (right), 
from which one can see clearly that the most massive haloes have 
nearly the same concentration in the two models, because the ef¬ 
fects of modified gravity are efficiently screened in these objects. 
It is well known from early studies that the halo concentration in 
ACDM simulations is given by a power-law function of mass (e.g., 
Bullock et al. 2001; Zhao et al. 2003; Neto et al. 2007; Duffy et al. 
2008; Maccio’ et al. 2008; Giocoli et al. 2010; Dooley et al. 2014), 
and our ACDM simulation shows the same result as illustrated by 
the red curves in Fig. 4 (neglecting the scatter at large halo masses, 
which is due to the small numbers of haloes there). Recent simu¬ 
lations and modelling have indicated that the mass dependence of 
the halo concentration can be more complicated and is not a simple 
powerlaw across the whole halo mass range (e.g., Prada et al. 2012; 
Sanchez-Conde & Prada 2014; Ludlow et al. 2014; Ng et al. 2014). 
However, our GR simulation has too small a dynamical range to be 
affected by this. 

In f(R) gravity, however, this is no longer true. Indeed, here 
we find a turning mass scale M,, below which the halo C-M 200 re¬ 


lation shows a clear deviation from a single power law and becomes 
higher than in GR. We have checked this discovery by running the 
Amiga Halo Finder (Knollmann & Knebe 2009, Ahf), which em¬ 
ployes a different method to measure halo concentrations, by using 
the relation between the maximum circular velocity and halo mass 
for NFW haloes, and found good agreement. We also make an addi¬ 
tional test by fitting the halo density profiles to the Einasto formula 
(Einasto 1965; Navarro et al. 2004), because it is known in ACDM 
that the shape of spherically averaged halo density profiles devi¬ 
ates systematically (though only slightly) from the two-parameter 
NFW formula and can be better described by the three-parameter 
Einasto formula (Gao et al. 2008). The Einasto fitting is less sensi¬ 
tive to the radius range used in the fitting, but in the test we only 
use radial bins outside the Power et al. (2003) convergence radius. 
Again, we have found very good agreement with Fig. 4. Finally, 
in Fig. 4 we have included all haloes, and in the last test we have 
also checked the results for relaxed haloes only, using the criteria 
proposed by Neto et al. (2007). We find that such a selection of re¬ 
laxed haloes makes very little difference in the concentration-mass 
relation, which agrees well with the findings of Gao et al. (2008). 
Since the main focus of this paper is the comparison between f(R) 
gravity and GR, we shall not show the plots from those tests. 

A possible reason for the difference in the concentration-mass 
relations of the two models studied here is the following: the turn¬ 
ing mass scale, AT*, which itself depends time, is roughly a thresh¬ 
old mass for the fifth force screening in f(R) gravity at each given 
time. Less massive haloes are unscreened and have deeper poten¬ 
tials than GR haloes with the same mass, which can make particles 
move towards the central regions and lead to higher concentrations. 
The increase of AT* with time reflects the simple fact that as time 
goes on more massive haloes become unscreened. 

Similar behaviour has also been found in other modified grav¬ 
ity theories. As an example, Barreira et al. (2014,?) find that, for 
models in which the strength of gravity increases rapidly in time, 
halos tend to be more concentrated (and vice versa). In chameleon- 
type theories, including f(R) gravity, the screening makes the sit¬ 
uation more complicated, but the general picture is that haloes tend 
to be more concentrated if the model has had an efficient screen¬ 
ing at early times (such as F6) because, at late times when screen¬ 
ing is ‘switched off', the potentials inside haloes deepen suddenly, 
and matter particles tend to fall towards the halo centre (Li & Zhao 
2010; Zhao et al. 2011a). Finally, in the phenomenological ReBEL 
model of Nusser et al. (2005), in which a scalar-mediated Yukawa- 
type fifth force helps in boosting the structure formation from early 
times, Hellwing et al. (2013) notice that the halo concentration is 
higher for all halo masses. These authors compare the kinetic and 
potential energies in their virialised haloes, and find that the ratio 
between the two is actually smaller than in ACDM haloes of same 
masses (cf. Fig. 11 in that paper). Even though the fifth force in the 
ReBEL model starts to effect from early times, the fact that it has a 
finite range (not longer than lft^Mpc in the models simulated by 
Hellwing et al. (2013)) means that the enhanced gravity could not 
affect regions beyond ~ l/i _1 Mpc: this is similar to the behaviour 
of the fifth force in F6 for our small haloes, which is possibly why 
the effect on the halo concentrations is also similar in the two cases. 

Another possible reason for the different c-M relations in F6 
and GR is the different halo formation histories in the two models. 
As is mentioned above, haloes which form at earlier times gener¬ 
ally have higher concentration because the mean matter density is 
higher when they collapse. Consider two (small) haloes of the same 
mass in GR and in F6: it is more likely that the latter has a larger 
fraction of its present-day mass assembled at later times, and thus 
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Figure 5. The halo formation time Zf as a function of M 2 oo, from our GR 
(red circles and curve) and F6 (blue triangles and curve) simulations. Error 
bars are I -rr standard deviations. 

its inner region is smaller, forms earlier and is more concentrated 
(in other words, a halo with mass Mi in F6 is likely to have a mass 
M 2 < Mi in GR and thus have a higher concentration than a GR 
halo of mass Mi). It would be useful to disentangle the two effects 
affecting halo concentrations, but this is difficult because a modi¬ 
fied gravitational force will always simultaneously affect both the 
halo accretion history and halo potential, except in cases where the 
screening is very strong inside haloes, such as in the cubic Galileon 
model (Barreira et al. 2014). We shall leave such a study for future 
work. 

We caution that the result for F6 may not quantitatively hold 
for HS f(R) models with other values of n or fno, or to other f(R) 
or chameleon models. The complicated physics that determines the 
concentration implies that the C-M 200 relation needs to be studied 
on a case-by-case basis in general. 

3.3 Halo formation histories 

The formation of dark matter haloes is a complicated process, in 
which frequent mergers and the accretion of smaller haloes hierar¬ 
chically lead to the formation of larger haloes. In this picture, large 
haloes form later when the environmental density is lower, and thus 
have lower concentrations than small haloes, as we have seen in the 
previous subsection. 

As the gravitational force is enhanced in f(R) gravity, is has 
been speculated that the matter clustering is stronger and as a result 
dark matter haloes form earlier in our f(R) model than in GR. For 
example, a previous study by Hellwing et al. (2010) found that in 
the ReBEL model, the Yukawa-type fifth force helps to form haloes 
at higher redshift than in the standard ACDM model, and there¬ 
fore can potentially move reionisation to earlier times as implied 
by CMB observations. Here, we want to study the halo formation 
times in our F6 simulations. 

In order to follow the growth of a halo with time, we start with 
the halo at the present time and identify its most massive progeni¬ 


tor from the previous snapshot. We repeat this procedure until the 
halo mass is too small to be resolved anymore, and define the halo 
formation time as the redshift, Zf, at which the most massive pro¬ 
genitor halo has assembled half of its mass at z = 0.0. This forma¬ 
tion time has been widely used in the literature (e.g., Lacey & Cole 
1993; Gao et al. 2004), although other definitions have also been 
used (e.g., Wechsler et al. 2002). 

In Fig. 5, we plot the halo formation time Zf as a function of 
M 200 ■ In both models, the results agree with the above hierarchi¬ 
cal picture that low-mass haloes form earlier. When comparing the 
two models, we can see that haloes more massive than 10 13 /i - 1 Mq 
form at nearly the same redshift in GR and F6, showing again that 
the chameleon mechanism works efficiently for these haloes to sup¬ 
press the effects of modified gravity. Less massive haloes, on the 
other hand, form slightly earlier in GR than in F6. This result seems 
to disagree with the general pattern found in Hellwing et al. (2010) 
in the ReBEL model, and seems in contrast to the naive expectation 
that the enhanced gravitational force in F6 boosts the hierarchical 
structure formation. 

To explain this behaviour, we need to again carefully examine 
the subtle differences between different models. In ReBEL, there is 
a Yukawa-type fifth force between particles, whose strength decays 
with distance but does not change in time. This implies that the fifth 
force starts to boost structure formation from early times, resulting 
in haloes forming earlier. In the case of F6, gravity is suppressed at 
redshift z > 1, and even at 2 < 1 it is only enhanced for smaller 
haloes. This means that: 

(i) the formation history of very massive haloes (e.g., M > 
10 13 /i _1 Mq) does not see the effect of an enhanced gravity, as we 
have seen above; 

(ii) less massive haloes evolve in a similar manner as in GR at 

z > 1, but grow more rapidly at 2 < 1, and as such they are more 
massive than their GR counterparts at present. As Zf is defined to 
be the time when a halo has gained half of its current mass (denoted 
by the halo would have a larger M 2 / 2 in F6 than in GR. But 

small haloes typically grow to M]y 2 at 2 > 1, before when there is 
little difference between GR and F6, and so it takes the halo longer 
to acquire a mass of M 1 / 2 in F6 than in GR, which means that the 
halo forms later in f(R) gravity. This, of course, is purely a result 
of the definition of Zf, and does not imply that matter clusters more 
slowly in F6. 

Therefore, like the concentration, the halo formation time also 
depends sensitively on the nature of gravity. Even for two models 
in both of which gravity is enhanced, the behaviour of c (M 2 00 ) or 
Zf (M 2 00 ) can be qualitatively different. For this same reason, the 
results for Zf for F6 can not be generalised to other variants of the 
HS f{R) model or other f(R) models without careful tests. 

3.4 Halo velocity dispersion profiles 

Before leaving this section, we study the velocity dispersion profile 
in our simulations, which is defined as 

& - ^) 2 > (14) 

P i£Ar 

in which i 6 A r means that particle i sits in a spherical shell from 
radius r — Ar/2 to r + Ar/2, and N p is the number of particles 
within this shell. Vi and Vh are the particle and host halo velocities 
respectively, and the latter is calculated as the average of the ve¬ 
locities of the 25% most bound particles in the host halo. The halo 
velocity dispersion is a more direct characterisation of the poten¬ 
tial inside a halo; it is determined by the dynamical (Schmidt 2010; 
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Figure 6. The stacked velocity dispersion profiles from our f(R) (filled triangles) and GR (filled circles) simulations, at three redshifts: z = 0 (left panel), 
z = 0.5 (middle panel) and z = 1 (right panel). We show the results for four bins of the host halo mass increasing from the bottom curves to the top curves, 
and indicated in the legend (we have divided the mass bin 10 12 ~ 10 13 /t —1 Mg into two sub-bins because this is the only bin which shows difference 
between F6 and GR). The solid curves simply connect the symbols, and the error bars show the l-er scatter around the mean. 


Zhao et al. 2011b) or effective (He et al. 2015) mass of a halo, and 
is enhanced by the modified force for unscreened haloes (Lee et al. 
2013; He et al. 2015). 

In Fig. 6 we show the velocity dispersion profiles measured 
from our simulations at z = 0.0 (left), 0.5 (middle) and 1.0 (right). 
Thanks to the chameleon screening, the difference between the two 
models for haloes more massive than ~ 10 13 /i _ 1 M© is almost un¬ 
detectable. Haloes in the mass range 10 12 — 10 can have 

significantly higher velocity dispersion in f(R) gravity than in GR, 
and the deviation increases with the distance from the halo centre, 
since the screening in f(R) gravity is relatively weak inside small 
halos, particularly in their outer regions in which matter density is 
low. We also notice that the enhancement of velocity dispersion is 
weaker at earlier times, due to stronger chameleon screening and 
less time for the fifth force to take effect. 

The result confirms that particles bound in unscreened haloes 
have higher kinetic energy to balance the extra potential produced 
by the fifth force. This implies that measurements of galaxy veloc¬ 
ity dispersions in galaxy groups, such as the Local Group, may not 
be able to give reliable estimates of the true masses of the systems. 
For example, to use such measurements to find the underlying mass 
requires a good understanding of the screening, which in turn re¬ 
quires an accurate knowledge of the true mass (as well as the envi¬ 
ronmental effects). Therefore, a trial-and-error procedure would be 
needed to improve the mass estimation iteratively from some initial 
guess, and each iteration needs to be calibrated by high-resolution 
simulations which take into account the full environmental effects 
and other complexities such as irregular shapes of haloes and dis¬ 
tributions of their massive satellites. 

On the other hand, if we indeed live in an unscreened region in 
f(R) gravity, but choose to interpret our measurements of galaxy 
velocity dispersions in the incorrect framework of GR, then the es¬ 
timated mass will be biased high compared with its true value. We 
will briefly mention one of its implications below. In any case, it is 
clear that f(R) gravity would make the already uncertain estimates 
of the Milky Way mass even more complicated. 


4 PROPERTIES OF SUBSTRUCTURES 

In the previous section we analysed the simulation results of vari¬ 
ous halo properties in our F6 and GR simulations. In this section, 
we turn our attention to the properties of subhaloes in these models. 

In hierarchical structure formation, halo merger events leave 
plenty of remnant structures that survive as subhaloes in the descen- 
dent haloes. As galaxies form inside haloes and migrate with them, 
subhaloes then exist as the host sites of satellite galaxies in galaxy 
groups and clusters. The properties of subhaloes and their evolu¬ 
tion history (i.e., the subhalo merger tree) provide the backbone for 
models of galaxy formation (see, e.g., Baugh 2006, for a review). 
The abundance and distribution of subhaloes also has important im¬ 
plications for the indirect detection of dark matter, for example by 
boosting the dark matter annihilation signal (e.g., Gao et al. 2012; 
Han et al. 2012). 

The fact that subhaloes form through hierachical mergers can 
also be utilised to identify them. Here we will use the tracking sub¬ 
halo finder Hierarchical Bound-Tracing (Han et al. 2012, HBT) to 
identify subhaloes. Starting from isolated haloes at an earlier snap¬ 
shot, HBT identifies their descendents at subsequent snapshots and 
keeps track of their growth. As soon as two haloes merge, HBT 
starts to track the self-bound part of the smaller progenitor as a 
subhalo in each subsequent snapshot. With a single walk through 
all the snapshots, all the subhaloes formed from halo mergers can 
be identified in this way. Such a unique tracking algorithm enables 
HBT to largely avoid the resolution problem suffered by configu¬ 
ration space subhalo finders (Muldrew et al. 2011; Han et al. 2012; 
Onions et al. 2012). By construction, HBT also produces clean and 
self-consistent merger trees that naturally avoid subtle defects such 
as missing links and central-satellite swaps common to many other 
tree builders (Srisawat et al. 2013; Avila et al. 2014). 


4.1 Subhalo mass functions 

Similar to haloes, the abundance of subhaloes can be described by a 
subhalo mass function (SHMF). The SHMF is known to depend on 
the size of their host haloes (Gao et al. 2004; van den Bosch et al. 
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Figure 7. The differential subhalo mass function dA/"/dM su b as a function of the subhalo mass M su b, at three redshifts - z = 0.0 (left panel), z = 0.5 
(middle panel) and 2 = 1.0 (right panel). The lower subpanels show the ratio between the F6 and GR results, and the binning scheme of host halo masses is 
indicated by different colours as shown in the legends (note that the highest mass bin does not exist in the right panels because at 2 : = 1 haloes more massive 
than 10 14 /i —1 M© does not exist in great numbers). The error bars (only shown in the lower subpanels for clarity) are l-cr standard deviations in each subhalo 
mass bin similarly as in the halo mass function plots. For guide eyes, we also plot the power-law fitting results of the subhalo mass functions using lines with 
the same colours (solid for GR and dashed for F6). 


2005) in ACDM simulations, but is close to a universal power law 
function of the subhalo mass, M su b, when normalised by the host 
halo mass, Mh os t • 

In Fig. 7 we plot the SHMF for three bins of host halo mass, 
10 12 ~ 10 13 /i _1 M 0 , 10 13 ~ 10 14 /i -1 M Q) 10 13 ~ 1O 14 /i _1 M 0 
(see the legends), at three redshifts, 2 = 0 (left), 0.5 (middle) and 1 
(right). For clarity, the results for the highest (lowest) mass bin are 
shifted upwards (downwards) by a decade. A quick visual inspec¬ 
tion of Fig. 7 indicates that the power-law relation holds true for 
the ACDM (circles and solid lines) and F6 (triangles and dashed 
lines) simulations as well, though the slope has a weak depen¬ 
dence on the host halo mass (lower for low-mass host haloes). To 
check this result, we tested HBT on a simulation using a different 
N-body code (described in Jing & Suto 2002) and found the same 
tendency. We also tested our simulations using the ROCKSTAR code 
(Behroozi et al. 2013) to identify subhaloes, but did not notice any 
dependence of this slope on the host halo mass. Therefore, we con¬ 
clude that this is likely due to the subhalo finding algorithm we 
use, which finds more massive and extended subhaloes than some 
other algorithms (Flan et al. 2012). We note that, even though the 
SF1MF from F1BT has a lower slope than the result from ROCK- 
STAR, it is consistently higher for the range of subhalo mass shown 
in Fig. 7. Because we are mainly interested in the relative differ¬ 
ences between models in this paper, we will leave a more detailed 
comparison of different algorithms to a future separate work and 
not show a plot for the comparison. 

From Fig. 7 we find that the difference between F6 and GR is 
smaller for more massive host haloes and at earlier times, because 
in both cases the chameleon screening is more efficient and effects 
of modified gravity more strongly suppressed. Differences between 
the two models also tend to be larger for small subhaloes, with F6 


predicting 20 ~ 50% more subhaloes with M su b between 10 11 and 
10 10 /i _ 1 Mq than GR in host haloes of mass 10 12 ~ 10 13 /i _ 1 M@. 
This implies that the enhanced gravity in the ,f(R) model studied 
here can help produce a substantially higher abundance of substruc¬ 
tures in Milky-Way-sized dark matter haloes. We will discuss the 
implication of this in the context of Milky Way satellite abundances 
below when discussing the subhalo velocity function. 

Note that an enhanced gravity will not only boost the cluster¬ 
ing of matter and formation of subhaloes, but can also increase the 
stripping of matter from subhaloes inside haloes (and thus decrease 
subhalo masses). Our results above suggest that the latter effect is 
subdominant. 

4.2 Subhalo spatial distributions 

Next we focus on the spatial distribution of subhaloes inside their 
host haloes. Naturally, one expects this distribution to depend on 
the nature of gravity,though this dependence can be weakened by 
the chameleon screening by the host haloes in f(R) gravity. 

Gao et al. (2004) showed that the spatial distribution of sub¬ 
haloes does not have a significant dependence on their host halo 
masses. In our ACDM simulations we have found the same result, 
as shown in Fig. 8, in which we plot the cumulative radial number 
distributions of subhaloes as circles for ACDM. We show in dif¬ 
ferent colours the results for three mass bins of host haloes, all at 
2 = 0, which agree well with each other. 

To see the effect of f(R) gravity, we also plot the correspond¬ 
ing results from the F6 simulation in Fig. 8 using triangles. There 
is very little difference from the GR results, possibly because of the 
efficient screening. Notice that here we have only shown results for 
host haloes more massive than 10 13 /i _1 Mq, in which the modified 
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Figure 9. The differential subhalo velocity functions from our f(R) (filled triangles) and GR (filled circles) simulations, at three redshifts: z = 0 (left panel), 
z = 0.5 (middle panel) and z = 1 (right panel). We show the results for three bins of the host halo mass increasing form the bottom curves to the top curves, 
and indicated in the legend; the solid curves simply connect the symbols. Note that the number of subhaloes in each bin is normalised by the total mass of the 
host halo. The results for the highest (lowest) host halo mass bin are artificially shifted upwards (downwards) by a decade to make the plot clearer. Error bars 
show the 1 -a scatter around the mean. 



Figure 8. The cumulative radial number distributions for subhalos, for three 
mass bins of their host haloes as indicated by the legends. The vertical axis 
is the fraction of subhaloes within r/r^oo, and we show the results for both 
F6 (triangles) and GR (circles) at z = 0. 

gravity effects are strongly suppressed as we have seen above. The 
results for smaller host haloes are not shown since they are noisier 
due to resolution limitations. 

4.3 Subhalo velocity function 

Subhaloes reside in the high-density environments within their host 
haloes, and experience constant tidal stripping, which strips mass 


from their outer parts. Their mass could change significantly dur¬ 
ing their evolution. Therefore, in the literature people often use the 
maximum circular velocity V max instead, because it depends pri¬ 
marily on the inner part of a (sub)halo. 

Following Gao et al. (2004), in Fig. 9 we plot the differential 
abundance of subhaloes as a function of V ma x, also known as the 
subhalo velocity function (SFIVF). The ACDM results in the range 
of 30 km/s < V max < 200 km/s are well described by a univer¬ 
sal power-law function, in agreement with findings in the literature 
(e.g., Gao et al. 2004; Dooley et al. 2014) (note that we have shifted 
the curves for different host halo mass bins for clarity), and drop off 
at small (large) Vmax is due to the resolution limit (finite box size). 

In f(R) gravity, the qualitative behaviour of the SHVF is sim¬ 
ilar, but the enhanced gravity leads to quantitative differences. For 
more massive host haloes, the difference is most significant at small 
Vmax, which correspond to smaller subhaloes that are less screened 
and therefore have formed in higher abundances; in contrast, larger 
subhaloes, with larger Vmax, are better screened and so their abun¬ 
dances do not change significantly from the GR predictions. For 
less massive host haloes, there is a noticeable boost in the subhalo 
abundance even for large subhalo Vnax, since the host haloes have 
become less screened since earlier times and substructures have 
more time to grow there. This dependence on host halo mass in 
principle implies a deviation from a universal SHVF, although the 
effect we see in Fig. 9 is fairly weak. 

The enhanced SHVF at V ma x > 30 km/s for host haloes with 
mass of ~ 10 12 hr 1 Mq seems to suggest that the missing mas¬ 
sive satellite problem of the Milky Way galaxy is worse in f(R) 
gravity, since in the latter the observed number of dwarf galaxies 
remains the same while the theoretically predicted number of mas¬ 
sive subhaloes is larger. Wang et al. (2012) argue that the missing 
satellite problem is not serious enough to motivate a revision to 
the ACDM paradigm, but what we saw above in Fig. 9 certainly 
seems to make ,f(R) gravity disfavoured. However, there are com¬ 
plicated issues which preclude a definite conclusion. For example, 
most measurements of the Milky Way halo actually predict its dy¬ 
namical mass, which can be 1/3 heavier than the true mass in f(R) 
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gravity - hence, with a given rotation speed of the Milky Way disk, 
the actual mass of the halo could be smaller than what we currently 
think. Also, galaxy formation can also be different in f(R) gravity, 
so that the way in which galaxies populate massive subhaloes might 
be different, making a direct comparison with ACDM even harder. 
We will leave detailed studies of these issues to future works. 


5 DISCUSSIONS AND CONCLUSIONS 

To briefly summarise, in this paper we have employed a very high- 
resolution simulation to study the properties of dark matter haloes 
and subhaloes of a f(R) gravity model. This model is a variant of 
that proposed by Hu & Sawicki (2007), with parameters n = 1 and 
| f hq | = 10 -6 . We argue that this is a borderline model that should 
be studied as a first step towards a more rigorous constraint on /.ro 
combining cosmological and astrophysical observations. We regard 
this as a realistic model which is not yet apparently ruled out by 
current data. 

The simulations we use in our analyses have 512 3 particles in 
a cubic box of size Lbox = 64/i _1 Mpc, with a background cos¬ 
mology chosen to be that of the best-fit WMAP9. This cosmology 
is more updated and more realistic than those of the previous f(R) 
simulations conducted by us (e.g., Li et al. 2013), and the resolu¬ 
tion here is also significantly higher, making it possible to study 
subhaloes in detail. Our halo catalogue is constructed using a stan¬ 
dard friend-of-friend linking method, and the subhaloes were found 
using the HBT algorithm of Han et al. (2012). 

Due to the efficient chameleon screening, this f(R) model 
shows small deviations from ACDM in general. For example, the 
halo mass function shows at most ~ 20% enhancement compared 
with the ACDM result between z = 0 and 2 = 1, with the deviation 
propagating to more massive haloes as time passes, in agreement 
with the semi-analytical predictions of Li & Efstathiou (2012). The 
dark matter distribution inside halos is almost identical in this f(R) 
model as in ACDM for haloes more massive than ~ 10 13 h^ 1 M@, 
again due to the chameleon screening; however, for smaller haloes, 
the screening is less efficient, which results in a deepening of the 
total potential and subsequently a steepening of the density profile. 
As a result, the halo concentration-mass relation is enhanced for 
such low-mass haloes and can no longer be described by a simple 
power law (as for ACDM). The stronger gravitational force in this 
f(R) model also enhances the growth of small haloes, but mainly 
at late times and as a result the halo formation time (i.e., the time 
by which a halo has gained half of its present-day mass) is actually 
later than in ACDM. We stress that these conclusions hold only 
for this specific f(R) model, and there is evidence suggesting that 
other models could behave qualitatively differently because of the 
complicated behaviour of gravity. We also notice enhanced halo 
velocity profiles in this f(R) model, confirming various previous 
work (e.g., Corbett Moran et al. 2015; Gronke et al. 2014; He et al. 
2015). 

The stronger gravity also helps to produce more substructures, 
mainly in host haloes less massive than ~ 10 13 h.~ 1 M@ because of 
the weaker screening therein, and for subhaloes less massive than 
10 12 h~ 1 M 0 . We find that Milky Way-sized haloes could host up to 
20 ~ 50% more subhaloes in the mass range 10 10 ~ 10 11 h~ 1 M & 
in the studied f(R) model than in ACDM. The subhalo mass func¬ 
tion can be fitted using a simple power law, as in ACDM, but with 
different parameters. We do not find a noticeable difference in the 
radial distribution of subhaloes inside their host haloes between the 
two models, though. The higher abundance of substructures is con¬ 


firmed in the subhalo velocity functions, which seems to make the 
missing satellite problem of the Milky Way worse. However, we 
stress that there are caveats in interpreting the result at its face, 
due to the further complexities in observationally determining halo 
mass in the context of modified gravity. 

Overall, we find that halo and subhalo properties of this bor¬ 
derline f(R) model are close to the ACDM predictions for mas¬ 
sive haloes, confirming previous results that this model is difficult 
to distinguish from ACDM using cosmological observations. How¬ 
ever, a substantial deviation might be found in less massive haloes 
such as that of our Milky Way, which is in agreement with the find¬ 
ings of previous low resolution simulations. This indicates that the 
dynamics of systems such as the Local Group can be sensitive to 
modifications of gravity of this kind and strength. This should be a 
focus of further studies in the future, following the recent progress 
in zoom simulations made by Corbett Moran et al. (2015). 

As mentioned above, this is a first step of a more detailed study 
of this borderline model, and here we have not touched the topic of 
astrophysical constraints, which is much more complicated. Studies 
of Jain et al. (2013) and Vikram et al. (2014) have demonstrated the 
potential of using astrophysical systems to improve the constraints 
on fn o. It would be useful to have a better understanding of the im¬ 
pact that environmental screening could have on those constraints. 
As in f(R) models the local behaviour of gravity usually depends 
on its environment at much larger scales, high-resolution or zoom 
simulations are important for calibrating the interpretation of as¬ 
trophysical observations. They are also important because they can 
provide more realistic quantifications of the environments for stel¬ 
lar evolution, which depends on the nature of gravitation sensitively 
(Davis et al. 2013). 

Obviously, improved constraints may or may not rule out this 
f(R) gravity model. However, with the progress in both numerical 
simulations and theoretical modelling, we are on a path towards 
better understanding. In such a sense, we are currently in the state 
of liminality 1 , and much effort is still in need to pass it. 
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the opening of a world of new knowledge. 



Haloes and subhaloes in f(R) gravity 13 


REFERENCES 

Arnold C., Puchwein E., Springel V., 2014a, Mon. Not. Roy. As- 
tron. Soc., 440, 833 

Arnold C., Puchwein E., Springel V., 2014b, (arXiv: 1411.2600) 
Avila S., Knebe A., Pearce F. R., Schneider A., Srisawat C., 
Thomas P. A., Behroozi P., Elahi P. J., Han J., Mao Y.-Y., Onions 
J., Rodriguez-Gomez V., Tweed D., 2014, Mon. Not. Roy. As- 
tron. Soc., 441, 3488 

Barreira A., Li B., Hellwing W. A., Baugh C. M., Pascoli S., 2014, 
J. Cosmo. Astropart. Phys., 1409, 031 
Barreira A., Li B., Hellwing W. A., Lombriser L., Baugh C. M., 
et al., 2014, J. Cosmo. Astropart. Phys., 1404, 029 
Baugh C. M., 2006, Reports on Progress in Physics, 69, 3101 
Behroozi P. S., Wechsler R. H., Wu H. Y., 2013, Astrophys. J., 
762, 109 

Bel J., Brax P., Marinoni C., Valageas P., 2014, (arXiv: 1406.3347) 
Bose S., Hellwing W. A., Li B., 2015. J. Cosmo. Astropart. Phys., 
in press 

Brax R, van de Brack C., Davis A.-C., Mota D. E, Shaw D. J., 
2007a, Phys. Rev., D76, 124034 

Brax R, van de Brack C., Davis A.-C., Mota D. F., Shaw D. J., 
2007b, Phys. Rev., D76, 085010 

Bullock J. S., Kolatt T. S., Sigad Y., Somerville R. S., Kravtsov 
A. V., et al., 2001, Mon. Not. Roy. Astron. Soc., 321, 559 
Cai Y.-C., Li B„ Cole S„ Frenk C. S„ Neyrinck M„ 2014, Mon. 
Not. Roy. Astron. Soc., 439, 2978 
Cai Y.-C., Padilla N„ Li B„ 2014, (arXiv:1410.1510) 

Carroll S. M., De Felice A., Duvvuri V., Easson D. A., Trodden 
M„ et al., 2005, Phys. Rev., D71, 063513 
Carroll S. M., Duvvuri V., Trodden M., Turner M. S., 2004, Phys. 
Rev., D70, 043528 

Clifton T., Ferreira P. G., Padilla A., Skordis C., 2013, Phys. Rep., 
513, 1 

Corbett Moran C., Teyssier R., Li B., 2015, Mon. Not. Roy. As¬ 
tron. Soc., 448, 307 

Davis A. C., Lint E. A., J. S„ Shaw D. J., 2013, Phys. Rev., D85, 
123006 

Davis M„ Efstathiou G.. Frenk C. S., White S. D. M„ 1985, ApJ, 
292, 371 

Dooley G. A., Griffen B. F., Zukin R, Ji A. P., Vogelsberger M., 
et al., 2014, Astrophys. J., 786, 50 
Duffy A. R„ Schaye J.. Kay S. T„ Dalla Vecchia C., 2008, Mon. 
Not. Roy. Astron. Soc., 390, L64 
Einasto J., 1965, Trudy Inst. Astrofiz. Alma-Ata, 51, 87 
Gao L„ Frenk C. S., Jenkins A., Springel V„ White S. D. M„ 2012, 
Mon. Not. Roy. Astron. Soc., 419, 1721 
Gao L„ Navarro J. F.. Cole S„ Frenk C., White S. D„ et al., 2008, 
Mon. Not. Roy. Astron. Soc., 387. 536 
Gao L., White S. D., Jenkins A., Stoehr F.. Springel V., 2004, 
Mon. Not. Roy. Astron. Soc., 355, 819 
Giocoli C., Tormen G., Sheth R. K., van den Bosch F. C., 2010, 
Mon. Not. Roy. Astron. Soc., 404. 502 
Gronke M., Llinares C., Mota D. F., Winther H. A., 2014, 
(arXiv: 1412.0066) 

Han J„ Frenk C. S., Eke V. R„ Gao L„ White S. D. M„ Boyarsky 
A., Malyshev D., Ruchayskiy O., 2012, Mon. Not. Roy. Astron. 
Soc., 427, 1651 

Han J.. Jing Y. P.. Wang H., Wang W., 2012, Mon. Not. Roy. As¬ 
tron. Soc., 427, 2437 

He J.-h., Hawken A. J., Li B., Guzzo L., 2015, 

(arXiv: 1501.00846) 


He J.-h., Li B., Hawken A. J., Granett B. R., 2014. Phys. Rev., 
D90, 103505 

He J.-h., Wang B„ 2013, Phys. Rev., D87, 023508 
Hellwing W. A., Cautun M., Knebe A., Juszkiewicz R., Knoll- 
mann S., 2013, J. Cosmo. Astropart. Phys., 10, 012 
Hellwing W. A., Knollmann S. R., Knebe A., 2010. Mon. Not. 
Roy. Astron. Soc., 408, 104 

Hellwing W. A.. Li B„ Frenk C. S., Cole S., 2013. Mon. Not. Roy. 
Astron. Soc., 435, 2806 

Hinshaw G., Larson D., Komatsu E., Spergel D. N., Bennett C. L., 
Dunkley J., Nolta M. R., Halpern M., Hill R. S., Odegard N., 
Page L., Smith K. M., Weiland J. L., Gold B., Jarosik N., Kogut 
A., Limon M., Meyer S. S.. Tucker G. S., Wollack E., Wright 
E. L., 2013, Astrophys. J. Suppl., 208, 19 
Hu W„ Sawicki I„ 2007, Phys. Rev., D76, 064004 
Jain B., Vikram V., Sakstein J., 2013, Phys. Rev., D79, 339 
Jenkins A.. Frenk C. S„ White S. D. M„ Colberg J. M„ Cole S„ 
Evrard A. E., Couchman H. M. R, Yoshida N., 2001, Mon. Not. 
Roy. Astron. Soc., 321, 372 

Jennings E., Baugh C. M., Li B., Zhao G. B., Koyama K., 2012, 
Monthly Notices of the Royal Astronomical Society, 425, 2128 
Jing Y„ Suto Y„ 2002, Astrophys. J„ 574, 538 
Joyce A.. Jain B., Khoury J.. Trodden M., 2014, arXiv: 1407.0059 
Khoury J.. Weltman A., 2004. Phys. Rev., D69, 044026 
Knollmann S., Knebe A., 2009, Astrophys. J. Suppl., 182, 608 
Lacey C., Cole S., 1993, Mon. Not. Roy. Astron. Soc., 262, 627 
Lee J., Zhao G.-B., Li B., Koyama K., 2013, Astrophys. J., 763, 
28 

Li B., Efstathiou G., 2012, Mon. Not. Roy. Astron. Soc., 421, 
1431 

Li B., Hellwing W. A., Koyama K., Zhao G., Jennings E.. Baugh 
C. M„ 2013, Mon. Not. Roy. Astron. Soc., 428, 743 
Li B., Lam T. Y., 2012, Mon. Not. Roy. Astron. Soc., 425, 730 
Li B., Zhao G., Koyama K., 2012, Mon. Not. Roy. Astron. Soc., 
421, 3481 

Li B., Zhao G., Teyssier R.. Koyama K., 2012. J. Cosmo. As¬ 
tropart. Phys., 01, 051 

Li B„ Zhao H„ 2010, Phys. Rev., D81, 104047 
Lombriser L., , 2014, Constraining chameleon models with cos¬ 
mology 

Lombriser L., Koyama K., Li B., 2014, J. Cosmo. Astropart. 
Phys., 1403, 021 

Lombriser L.. Li B., Koyama K.. Zhao G.-B., 2013, Phys. Rev., 
D87, 123511 

Lombriser L., Simpson F., Mead A., 2015, (arXiv:1501.04961) 
Ludlow A. D., Navarro J. F.. Angulo R. E., Boylan-Kolchin M., 
Springel V., et al., 2014, Mon. Not. Roy. Astron. Soc., 441, 378 
Maccio’ A. V., Dutton A. A., Bosch F. C. d., 2008, Mon. Not. Roy. 
Astron. Soc., 391, 1940 

Muldrew S. I., Pearce F. R., Power C., 2011, Mon. Not. Roy. As¬ 
tron. Soc., 410, 2617 

Navarro J. F. Frenk C. S., White S. D., 1996, Astrophys. J., 462, 
563 

Navarro J. F.. Frenk C. S., White S. D., 1997, Astrophys. J., 490, 
493 

Navarro J. F., Hayashi E., Power C., Jenkins A., Frenk C. S., et al., 
2004, Mon. Not. Roy. Astron. Soc., 349, 1039 
Neto A. F., Gao L.. Bett P., Cole S., Navarro J. F., et al., 2007, 
Mon. Not. Roy. Astron. Soc., 381, 1450 
Ng K. C. Y., Laha R., Campbell S., Horiuchi S., Dasgupta B., 
Murase K., Beacom J. F., 2014, Phys. Rev. D, 89, 083001 
Nusser A., Gubser S. S., Peebles R, 2005, Phys. Rev., D71, 



14 Shi et. al. 


083505 

Onions J., Knebe A., Pearce F. R., Muldrew S. I., Lux H., Knoll- 
mann S. R., Ascasibar Y., Behroozi P., Elahi P., Han J., Ma- 
ciejewski M., Merchan M. E., Neyrinck M.. Ruiz A. N., Sgro 
M. A.. Springel V.. Tweed D., 2012, Mon. Not. Roy. Astron. 
Soc., 423, 1200 

Power C., Navarro J., Jenkins A., Frenk C., White S. D., et al., 
2003, Mon. Not. Roy. Astron. Soc., 338, 14 
Prada F., Klypin A. A., Cuesta A., Betancort-Rijo J. E., Primack 
J., 2012, Mon. Not. Roy. Astron. Soc., 423, 3018 
Prunet S., Pichon C., Aubert D., Pogosyan D., Teyssier R.. Gott- 
loeber S., 2008, Astrophys. J. Suppl., 178, 179 
Reed D., Bower R., Frenk C., Jenkins A., Theuns T., 2007, Mon. 
Not. Roy. Astron. Soc., 374, 2 

Sanchez-Conde M. A., Prada F., 2014, Mon.Not.Roy.Astron.Soc., 
442, 2271 

Sawala T., Frenk C. S., Crain R. A., Jenkins A.. Schaye J., Theuns 
T., Zavala J., 2013, Mon. Not. Roy. Astron. Soc., 431, 1366 
Schaller M., Frenk C. S.. Bower R. G., Theuns T., Jenkins A., 
et al., 2014, (arXiv: 1409.8617) 

Schmidt F„ 2010, Phys. Rev., D81, 103002 
Sheth R. K„ Tormen G., 1999, MNRAS, 308, 119 
Srisawat C., Knebe A., Pearce F. R., Schneider A., Thomas P. A., 
Behroozi P, Dolag K., Elahi P. J., Han J., Helly J., Jing Y., Jung 
I., Lee J., Mao Y.-Y., Onions J., Rodriguez-Gomez V., Tweed D., 
Yi S. K„ 2013, Mon. Not. Roy. Astron. Soc., 436, 150 
Taruya A., Nishimichi T., Bernardeau F., Hiramatsu T., Koyama 
K„ 2014, Phys. Rev., D90, 123515 
Teyssier R., 2002, Astron. Astrophys., 385, 337 
van den Bosch F. C., Tormen G., Giocoli C., 2005, Mon. Not. Roy. 
Astron. Soc., 359, 1029 

Vikram V., Sakstein J., Davis C., Neil A., 2014, 

(arXiv: 1407.6044) 

Wang J., Frenk C. S., Navarro J. E, Gao L., Sawala T., 2012, Mon. 
Not. Roy. Astron. Soc., 424, 2715 
Wechsler R. H., Bullock J. S., Primack J. R., Kravtsov A. V., 
Dekel A., 2002, Astrophys. J., 568, 52 
Weinberg D. H., Mortonson M. J., Eisenstein D. J., Hirata C., 
Riess A. G., Rozo E., 2013, Phys. Rep., 530, 87 
White M., 2001, Astron. Astrophys., 367, 27 
Zhao D. H., Jing Y. P. Mo H. J., Borner G., 2003, Astrophys. J., 
597, L9 

Zhao G., 2014. Astrophys. J. Suppl., 211, 23 
Zhao G„ Li B„ Koyama K„ 2011a, Phys. Rev., D83, 044007 
Zhao G.. Li B„ Koyama K„ 201 lb, Phys. Rev. Lett., 107, 071303 
Zivick P., Sutter P, Wandelt B. D., Li B., Lam T. Y., 2014, 
(arXiv:1411.5694) 



